function data2fftMaySTHarmonic
format long;

Mayx50=load('Mayx50.out');
Mayx100=load('Mayx100.out');
Mayx150=load('Mayx150.out');
Mayx200=load('Mayx200.out');
Mayx250=load('Mayx250.out');
Mayx300=load('Mayx300.out');
Mayx350=load('Mayx350.out');
Mayx400=load('Mayx400.out');
Mayx450=load('Mayx450.out');
Mayx500=load('Mayx500.out');
Mayx550=load('Mayx550.out');
Mayx600=load('Mayx600.out');
Mayx650=load('Mayx650.out');

datax50=[Mayx50(:,1) Mayx50(:,2)];
datax100=[Mayx100(:,1) Mayx100(:,2)];
datax150=[Mayx150(:,1) Mayx150(:,2)];
datax200=[Mayx200(:,1) Mayx200(:,2)];
datax250=[Mayx250(:,1) Mayx250(:,2)];
datax300=[Mayx300(:,1) Mayx300(:,2)];
datax350=[Mayx350(:,1) Mayx350(:,2)];
datax400=[Mayx400(:,1) Mayx400(:,2)];
datax450=[Mayx450(:,1) Mayx450(:,2)];
datax500=[Mayx500(:,1) Mayx500(:,2)];
datax550=[Mayx550(:,1) Mayx550(:,2)];
datax600=[Mayx600(:,1) Mayx600(:,2)];
datax650=[Mayx650(:,1) Mayx650(:,2)];

NNN=length(Mayx50);

fftx50=fft(datax50(:,2))/NNN;
fftx100=fft(datax100(:,2))/NNN;
fftx150=fft(datax150(:,2))/NNN;
fftx200=fft(datax200(:,2))/NNN;
fftx250=fft(datax250(:,2))/NNN;
fftx300=fft(datax300(:,2))/NNN;
fftx350=fft(datax350(:,2))/NNN;
fftx400=fft(datax400(:,2))/NNN;
fftx450=fft(datax450(:,2))/NNN;
fftx500=fft(datax500(:,2))/NNN;
fftx550=fft(datax550(:,2))/NNN;
fftx600=fft(datax600(:,2))/NNN;
fftx650=fft(datax650(:,2))/NNN;

% Write the initial conditions to a file.
% These initial conditions are needed by the Fortran solvers.

ZeroModeLoc=181;
deltak=1/(datax50(NNN,1)-datax50(1,1));
kvec=transpose(0:deltak:(NNN-1)*deltak);
k0=1.52806403275206526
epsilon=2*k0*abs(fftx50(ZeroModeLoc))

ICs=[transpose(kvec); transpose(real(fftx50)); transpose(imag(fftx50))];
fileID = fopen('ICST.out','w');
fprintf(fileID,'%20.16f %20.16f %20.16f\n',ICs);
fclose(fileID);

% Compute M norms and exponential fit.

x50Mnorm=0;
x100Mnorm=0;
x150Mnorm=0;
x200Mnorm=0;
x250Mnorm=0;
x300Mnorm=0;
x350Mnorm=0;
x400Mnorm=0;
x450Mnorm=0;
x500Mnorm=0;
x550Mnorm=0;
x600Mnorm=0;
x650Mnorm=0;

for j=136:226
    x50Mnorm=x50Mnorm+abs(fftx50(j))^2;
    x100Mnorm=x100Mnorm+abs(fftx100(j))^2;
    x150Mnorm=x150Mnorm+abs(fftx150(j))^2;
    x200Mnorm=x200Mnorm+abs(fftx200(j))^2;
    x250Mnorm=x250Mnorm+abs(fftx250(j))^2;
    x300Mnorm=x300Mnorm+abs(fftx300(j))^2;
    x350Mnorm=x350Mnorm+abs(fftx350(j))^2;
    x400Mnorm=x400Mnorm+abs(fftx400(j))^2;
    x450Mnorm=x450Mnorm+abs(fftx450(j))^2;
    x500Mnorm=x500Mnorm+abs(fftx500(j))^2;
    x550Mnorm=x550Mnorm+abs(fftx550(j))^2;
    x600Mnorm=x600Mnorm+abs(fftx600(j))^2;
    x650Mnorm=x650Mnorm+abs(fftx650(j))^2;
end
 
positions=epsilon^2*k0*[0,50,100,150,200,250,300,350,400,450,500,550,600]';
xdata=epsilon^2*k0*linspace(0,650,100)';
Mnorms=(k0/epsilon)^2*[x50Mnorm,x100Mnorm,x150Mnorm,x200Mnorm,x250Mnorm,x300Mnorm,x350Mnorm,x400Mnorm,x450Mnorm,x500Mnorm,x550Mnorm,x600Mnorm,x650Mnorm]';
logMnorms=log(Mnorms);
coeffs=polyfit(positions,logMnorms,1);
Mpred=exp(polyval(coeffs,xdata));
delta=-coeffs(1)/2
L2Scaling=exp(2*delta.*(positions));
ScaledMnorms=Mnorms.*L2Scaling;

figure(3)
subplot(2,1,1)
hold off
scatter(positions,Mnorms,'filled','k')
hold on
plot(xdata,Mpred,'r','LineWidth',2)
text(0.05,0.08,strcat('Computed \delta=',num2str(delta)),'FontSize',13,'FontName','Times')
ylabel('$\mathcal{M}$','Interpreter','Latex')
legend('data','best fit exp','Location','NorthEast')
title('$\mathcal{M}$','Interpreter','Latex')
axis tight

subplot(2,1,2)
hold off
scatter(positions,ScaledMnorms,'filled','k')
hold on
plot([0 0.33],[(k0/epsilon)^2*x50Mnorm (k0/epsilon)^2*x50Mnorm],'r','LineWidth',2)
text(0.05,0.47,strcat('Using \delta=',num2str(delta)),'FontSize',13,'FontName','Times')
xlabel('$\chi$','Interpreter','Latex')
ylabel('$\mathcal{M}$','Interpreter','Latex')
leg1=legend('scaled data','$\mathcal{M}(0)$','Location','NorthEast');
set(leg1,'Interpreter','Latex');
title('Scaled $\mathcal{M}$','Interpreter','Latex')
axis tight
